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Abstract. Based on reasonable testing model problems, we study the 
preservation by symplectic Runge-Kutta method (SRK) and symplectic 
partitioned Runge-Kutta method (SPRK) of structures for fixed points 
of linear Hamiltonian systems. The structure-preservation region pro- 
, vides a practical criterion for choosing step-size in symplectic computa- 

tion. Examples are given to justify the investigation. 

1. Introduction and preliminaries 
Consider the n-degree-of-freedom (d.o.f) Hamiltonian system 

I 



-I 



> • (1.1) z = JV.ff(z), J 

where H is a smooth scalar function of the extended phase space variables 

■ z 6 R , denoting the Hamiltonian, and J is the Poisson matrix with / the 
\ n x n identity matrix. By introducing the canonically conjugate variables, 

O ' z= (QjP)i t ne above system can be rewritten as 

OO ' 

(1.2) q = dH/dp, p = -dH/dq, 

where q £ M. n represents the configuration coordinates of the system and 
their canonically conjugate momenta p £ W 1 represents the impetus gained 

■ by movement. As is well-known, Hamiltonian systems are introduced as a 
type of system for which the existence of conservative quantities are auto- 
matic. System (|1 .2|) possesses two remarkable properties: 

(1) the solutions preserve the Hamiltonian, i.e., 

(2) the corresponding flow is symplectic, i.e., 
(1.4) j t [dpAdq] = 0. 
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In the last two decades, enormous attention has been paid to numerical 
methods which preserve the symplecticity, namely, symplectic integrators 
for Hamiltonian systems; we refer to the monographs Hairer et al. [1\ and 
Sanz-Serna&Calvo [8] for details and related literature. Theoretical analy- 
sis together with numerous numerical experiments has shown that symplec- 
tic integrator not only produces improved qualitative numerical behaviors, 
but also allows for a more accurate long-time scale computation than with 
general-purpose methods. In the symplectic integration study, a widely rec- 
ognized fact is that the symplecticity of a numerical integrator has little to do 
with its step-size. Particularly, for SRK and SPRK methods, their symplec- 
ticities are only related to the coefficients (see Section 2 below). Therefore, 
in practical computations, one usually resorts to the classical stability anal- 
yses to find a suitable range for choosing numerical step-sizes. However, in 
a recent paper it is shown that in some cases even the step-size of the 
symplectic Euler method satisfies the classical linear stability requirements, 
one can still get periodic- two numerical solutions, or even chaotic solutions. 
That means, we need to require more stringent conditions on step-sizes in 
addition to the classical stability considerations in simulations of Hamilton- 
ian flows, even for symplectic integrators. In the present paper, we make a 
first step towards such investigation by studying the influence induced by 
the numerical discretization on the equilibria structure of the underlying 
Hamiltonian system. It is recalled that for a general ODE of the form 



it may admit the presence of equilibrium point, namely, z G R m such that 
/(z) = 0, and the eigenvalues of the corresponding stability matrix V z /(z) 
determine the type of the equilibrium point and its stability properties. 

In the sequel, we are mainly concerned with the Runge-Kutta (RK) meth- 
ods and partitioned Runge-Kutta (PRK) methods. Henceforth, we cus- 
tomarily refer to an s-stage RK method by the triple 1Z S = (A,b,c), with 
A = (<Hj)fj =1 , b = (bi)f =1 and c = (cj)f =1 being, respectively, the coeffi- 
cient matrix, weights and abscissae, and an s-stage PRK method by the 
pair n { s ] -nf \ Next, we would like to review some of the classical linear 
stability concepts and by tracing the origins we can thus draw forth our mo- 
tivations for the current work. The probably most well known A-stability is 
introduced by Dahlquist in 1960's (see, e.g., [2j). Applying TZ S to the famous 
Dahlquist test equation 



z = 



/(z), zeK m , / : R m h-> R ! 



711 




and z = Xh 



with R(z) the stability function of 1Z S (see, Chapter IV. 3, [2]). It is noted 
that the solution to (|1.5j) asymptotically decays to zero as t — ► oo, and 
in order for the numerical scheme (|1.6p to yield such qualitative behavior 
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without any restriction on the step size h, we naturally require that 
(1.7) \R(z)\<l, forany/i>0. 

Methods satisfying (II. 7h are called A-stable, and this concept has been play- 
ing an indispensable role in the numerical field. Apparently, one can derive 
the same conclusion (II. 7p for 1Z S when applying it to the following equation 

(1-8) y' = Ay, 

where A £ C is the complex conjugate to A in equation (jl.5j> . If we set 
A = a + i/3 with a, [3 G R and a < 0, it is easy to see that equations (II. 5p 
and (jl.8l) are equivalent to the following system of ODE, 



(1.9) 



x = ax — /3y, 
y = j3x + ay. 



System (|1.9p has an equilibrium point (0, 0) and its corresponding stability 
matrix is given by 



which has two eigenvalues A 



(1.10) 
with 



Xl+l 

yi+i 



1,2 



Q 



a -13 
(3 a 

a ± i/3. Now, we apply 1Z S to (|1.9p and get 



R(z) 
R(z) 



Q 



-i 



z = Xh, z = Xh, X = a + i/3 and Q 



V2 



xi 

yi 

i 
1 



Introducing the forward difference operators 



St xi 



xi+i - XI 

h • 



scheme (ll.lOp can be reformulated into 
(1.11) t 



xi 
JJi 



Q 



R{z)-1 
h 







R(z)-1 



yi+i - yi 
h • 



Xl 

JJI 



which is the discrete dynamical system approximating (|1.9p . Obviously, 
(0, 0) is the fixed point (equilibrium point) for (|l.lip and the correspond- 
ing stability matrix is given by the coefficient matrix in (II. lip , which is 
seen to possess two eigenvalues, A^i = (R(z) — l)/h,\h t 2 = Xh,l- Now, it 
is readily seen that an A-stable RK method preserves the equilibria struc- 
ture unconditionally after discretization. In the reverse, if we want an RK 
method to give rise to certain preservation of the equilibria structure of the 
underlying ODE, we are naturally led to the condition (|1.7p . From above 
analyses, we can further easily deduce that the so-called stability region for 
an RK method (see, Chapter IV. 3, [2]) is indeed the set of those step sizes 
with which the RK method can preserve the equilibria structure of the test 
equation (\l.9h . 
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Next, we briefly mention one more example to illustrate the close con- 
nection between the so-called P-stability and the preservation of equilibria 
structure. In the theory of orbital mechanics, many problems are formulated 
as some special second order ODE 

d 2 u 

(1.12) — = f(t,u), u(t ) = u , u(t ) = u' , 

which is often known in advance to have periodic solution. In [3], Lambert 
and Watson introduced the concept of P-stability for numerical methods to 
solve equation (|1.12[) . It is remarked that in the P-stability is defined 
for multistep methods, while it is equally applicable to RK method. The 
test model problem given in [3] is 

(1.13) u = -A 2 u, AeR. 

Set u = v, then equation f)l . 13|) can be reformulated into a first order system 



(1.14) 



u = v, 
v = —X 2 u. 



In like manner as above for A-stability, one can show that the P-stable RK 
methods for (|1.13|) preserve the equilibria structure of (|1.14p unconditionally 
and vice versa. However, we will not explore more details on this aspect and 
now turn our study to Hamiltonian system (jl.ip . 

The equilibrium points of system (jl.ip are those (p, q) G W 1 x R n such 
that VH(p, q) = and the corresponding stability matrices are given by 

d(-dH/dq,dH/dp) ,_ 

JsM ■= ^ (p,«). 

To ease our study, we start with linear Hamiltonian systems, and then pro- 
pose several ways to extend our analyses to nonlinear case. We first desig- 
nate some model problems for our investigation. Apparently, these model 
problems should feature the Hamiltonian systems. In [7], it is shown that 
an n d.o.f linear Hamiltonian system can be canonically transformed into a 
Hamiltonian system consisting of n 1-d.o.f subsystems and these subsystems 
assume the following standard forms, 



(1.15) 



p = -P\ 

q= p, (3>0 



or 

J p = -/3p, 



(1.16) 



/3q, P > 0. 



Therefore, it is natural for us to take systems (I1.15P and (11.161) as our testing 
model problems and which will be used throughout. Clearly, (0, 0) is their 
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equilibrium point, and for system (j!.15[) 

-fi 2 

1 

with corresponding eigenvalues Ai,2 = which are of elliptic type, and we 
call system (|1 . 15[) have elliptic equilibria structure; while for system ()1.16|) 

-13 


with the eigenvalues = ±/3, which are of hyperbolic type, and we call 
system f 1 1 . 1 6 1) having hyperbolic equilibria structure. Next, we shall inves- 
tigate under what conditions, the equilibria structure can be inherited by 
the corresponding numerical schemes. And what would be caused by the 
numerical discretization if the structure-preservation conditions are ruined. 

In the section followed, we present some necessary and sufficient condi- 
tions for an SRK/SPRK method and its corresponding symmetric compo- 
sition method to preserve the equilibria structures of the testing problems. 
The fundamental idea is that we regard the numerical scheme as a discrete 
dynamic which is an approximation to the Hamiltonian system and there- 
fore should have the same equilibria structure as its continuous counterpart. 
Results show that in some cases, SRK method can preserve the equilibria 
structure unconditionally, but without exception, PRK method always gives 
a structure-preservation region on the positive real line for the numerical 
step-size. Furthermore, it is found that if the order of accuracy of a SPRK 
method is increased by symmetric composition, the structure-preservation 
region will be decreased accordingly. 

In section 3, we extend the analyses to nonlinear Hamiltonian systems. 
Besides, we give an example which justifies the necessities of such equilibria 
structure preservation study. 

2. Equilibria structure preservation of SRK/SPRK method 

In this paper, we confine our study to symplectic Runge-Kutta (SRK) 
methods and symplectic partitioned Runge-Kutta (SPRK) methods. The 
symplecticity conditions for 1Z S and IZ^p — Tig are, respectively, given by 

(see mm) 

(2.1) BA + A T B - bb T =0, B = diag[6]; 

(2.2) B ^AW+A^ T B^-b^b^ T =0, flW = diag[6«] (i = 1,2), 

(2.3) 6«=6< 2 ). 

2.1. Preserving the elliptic equilibria structure. We first apply 1Z S to 
system (j!.15j) and get the following scheme 



(2.4) 



P = pie s - hp 2 AQ, Pl+1 = Pl - hp 2 b T Q, 
Q = Qie s + hAP, qi +1 =q t + hb T P, 
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where P = [Pi,P 2 , . . . , P S ] T and Q = [Q u Q 2) . . . , Q s f are internal stage 
values, e s = [1,...,1] T and p\ p(lh),qi « q(lh). It can be computed 
straightforwardly that 

f P = (I + / l 2 /3 2 ^ 2 )~ 1 (e^ " ty? 2 ^), 
\ Q = (/ + ^ 2 A 2 r 1 ( es g i + Me 8 |),) 1 

and thus scheme (|2.4p is equivalently reformulated as 

?2 n 4 „, _ o2 



(2.5) 



= —h(3 2 DAe s pi - p 2 De s qi, 
S^qi = De s pi - hf3 2 DAe s q h 



where D = b T (I + h 2 f3 2 A 2 ) 1 and the forward difference operators are de- 
fined as 

Pl+i ~ Pi r + _ Ql+l ~ 1i 



h ' 1 h 
(|2.5p is a discrete dynamical system which approximates (jl,15p and obvi- 
ously, (0,0) is its fixed point (equilibrium point). The stability matrix of 
(|23D is 

-h(5 2 DAe s -P 2 De s 
De s -hj3 2 DAe s 



and its eigenvalues are 



A 



1,2 



-h(5 2 DAe s ±i(/3De 9 ), 



which are of elliptic type. Hence, the SRK method preserves the elliptic 
equilibria structure of system f)l . 15|) unconditionally, but here we note that 
there is a little shift of Ai^ from Ai ; 2 due to the numerical discretization. 

Next, we apply an s-stage SPRK method TZ^P — to system f)1.15j) and 
get the following scheme 



(2.6) 



P = Pl e s -hp 2 A^Q, 
Q = gi e s + hAWp, 



pi - h(3 2 b^ T Q, 
+ hb^ T P. 



Pi+i 

m+i = qi 

The discrete dynamical system equivalent to (12. 6ft is given by 



(2.7) 



6+ Pl = -hp 2 D^A^e sPl - l3 2 DWe sgi , 
S+ qi = D^e sPl - h[3 2 D^A^e s qu 



where 

(2.8) £>M = b^ T (I+h 2 p 2 A^A^y\ 
The stability matrix for ()2.7[) is 



. h p2 D (l) A (2) es 



L>(2) =b^ T (I+h 2 (3 2 A^A^)-\ 



-l3 2 D^e s 

. h p2 D (2) A (l) es 
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and its eigenvalues are 



Xi,2=U -hp 2 [D^A^e s +D^A^e s 



± yh 2 /P[DW>A( 2 )e a - D( 2 )iWe s ] 2 - A/3 2 [(DWe s )(D^e s )} 

Therefore, if the equilibria point (0, 0) for (|2.7p is of elliptic type we need to 
require 

(2.9) h 2 (5 2 [D^A^e s - D^A^e s ] 2 - 4[(D^e s )(D^e s )] < 0. 

Next, we introduce 

_ det(J + z 2 AW(AW - e a bW T )) 

, s ai ''~ det( I + z 2 A( 2 ) AW) 

2-10 V ' _ 

_ det(I + z 2 AW(AM - e s bM T )) 
a2 ''~ det(I + z*A0-)A<?)) 

where z = f3h, which shall also be used throughout, then condition (|2.9p can 
be reformulated as 

(2.11) \a 1 + a 2 \<2. 

In fact, from (I2.7p . one can get 

Pl+l =(1 - z 2 D^A^e s )pi - 0zD^e s q h 

qi + i =hD^e sPl + (1 - z 2 D^A^e s ) qi . 

Since the method is symplectic, it has 

dpi +1 A dqi +1 = dpi A dq h 

Substituting the relations (|2.12p into the above equality, then through straight- 
forward calculations, we obtain 

(D^e s )(D^e s ) = [D& A® e. + D® A® e 8 ] - z 2 (D^ A™ 'e s )(D^ A® <e s ), 

which is then substituted into (|2.9p to yield 

(2.13) |1 - z 2 D^A^e s + 1 - z 2 D^A^e s \ < 2. 

Next, by observing that 

1 _ z 2 D^A^e s =1 - z 2 b^ T (I + h 2 p 2 A^A^)- 1 A^e s = a x , 

1 _ z 2 D^A^e s =1 - z 2 b^ T (I + h 2 p 2 A^A^)- 1 A^e s = a 2 , 

we finally arrive at the equivalence of (|2.9p and (|2.1ip 
In summary, we have 
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Proposition 2.1. For system U.15\) . an SRK method 1Z S preserves its el- 
liptic equilibria structure unconditionally, whereas for an SPRK method, 
IZs — TZs , its structure-preservation region is given by 

(2.14) {h > : |ai +a 2 \ < 2}, 
where a% and a 2 are defined by \2.1(A) . 

Consider the symplectic Euler method, where = 0, = 1 and 
£>(!) = fr( 2 ) = 1. By ([230]) . it is computed that d = 1 - z 2 ,a 2 = 1, and 
hence for symplectic Euler method to preserve the elliptic structure of sys- 
tem (fTTT5l) we must have from (pTTIj) that |2 - z 2 \ < 2, i.e., < (3h < 2. 

Using the W-transformation due to Hairer and Wanner ( see Chapter IV. 5, 
[2]), we have 

(2.15) W~ X A^W = X A{1) , W~ l A^W = X A(2 ) , 

where W is the generalized Vandermonde matrix, and X A (i)(i = 1,2) are 
matrices possessing some standard form. In particular, for a class of impor- 
tant symplectic method, Lobatto IIIA-Lobatto IIIB pair (see [9~1 I10j). 



(2.16) 



X 



AW 



-a 





-6-2 

o 

6-i 



\ 







(2.17) 



X 



A(2) 



-6 
o 



V 



o 







/ 



where £ fc = l/2ViF~ 
using (|2.15p . Oi(i = 1 



(2.18) 



(2.19) 



d 



« 2 



- 1, = 1, 2, . . . , s. In term of W-transformation, i.e., 
2) in (|2.10p can be reformulated as 

_det(/ + z 2 X 4(2) (X 4(1) 



eV T )) 



det(I + z 2 X^2X j4 (i) 
_det(J + z 2 X A(1) (X A( 2 



eV T )) 



det(I + z 2 X A iX A(2) ) 

where e 1 = [1,0,..., 0, 1] T . For Lobatto IIIA-Lobatto IIIB pair, with X A a) (i - 
1,2) given in (12. 16ft and (12.17p . it is easily verified that d = a 2 and there- 
fore, for Lobatto IIIA-Lobatto IIIB method to preserve the elliptic structure 
of system (I1.15P we must have from (I2.14|) that |ai| < 1. In the follow- 
ing, we list some of the computed results for Lobatto IIIA-Lobatto IIIB 
method in term of its stage s; see Table [TJ Moreover, for s = 5, the 
structure-preservation region is < z < 3.140328, and for s = 10, it is 
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Table 1. Structure-preservation region for Lobatto IIIA- 
Lobatto IIIB pair 



s 




Z 


2 




< z < 2 


3 


x 24^ ^48 
1+ J- Z 2 


< z < \/8(~ 2.828) 


4 


24 

15 ^ ^ 900 3000 - 


< z < \/42 - 6y/29(^ 3.1127) 





< z < 3.141590. We are naturally led to the conjecture that as s — > oo, 
the elliptic structure-preservation region for an s-stage Lobatto IIIA-Lobatto 
IIIB method is given by < z < it. 

Next, we consider the composition of a given basic one-step method with 
different step sizes, 

(2.20) V h = $ ysh o • • • o $ 7lh , 

where $h 1S the basic SRK method and 7j/t(i = 1, . . . , s) are the composition 
step length (cf. [6].[Hj)- The composition is required to be symmetric, i.e., 
ji = 7<j_|_i_j(i = 1,2,..., [s/2]), and hence ^>h is still an SRK method. By 
Proposition ^. 11 we know that the symmetric composition of an SRK method 
still preserves the elliptic structure of system (|1.15p unconditionally. Next, 
consider the structure-preservation region of a symmetric composition of 
a SPRK method, which is known to be again a SPRK method. To ease 
our study, we only take the 2-stage Lobatto IIIA-Lobatto IIIB pair as an 
example, whose 4th-order symmetric composition is given in [5], and for 
which we compute 

oi = 1 " \z 2 + \llll^>li + ^l2)z A - ^7?72(7i + 12? 'z 6 , 
where 71 = 1/(2 - 2 1 / 3 ) and 72 = 2 1 / 3 /(2 - 2 1/3 ). By |oi| < 1, we get 

< z < v / 2l8(w 1.5748). 

Similarly, we further computed that the structure-preservation region for the 
corresponding 6th-order composition method is given by < z < 1.1034. 
Noting by Table [U the structure-preservation region for the underlying basic 
method is < z < 2. Base on this example, we conjecture that as the order 
of a SPRK method is increased by symmetric composition, its structure- 
preservation region will be decreased accordingly. A stringent proof of such 
conjecture is fraught with difficulties and is beyond the scope of this paper. 

2.2. Preserving the hyperbolic equilibria structure. The scheme of 
an SRK method 1Z S for system (|1 . 16|) is read as 

\ r+ R(Z)-1 
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where R(z) = det(/ + zA)/ det(I — zA) is the stability matrix for 1Z S . It is 
readily observed that for the discrete dynamical system (|2.2ip to preserve 
the equilibria structure of system (11 . 16f) . we must have 

(2.22) R(-z) - 1 < 0, and R{z) - 1 > 0, 

i.e., the hyperbolic structure-preservation region for an SRK method 1Z S is 

det(I-zA) det(I + zA) 

Since we always have R(—z) < l(z = (3h > 0) for an A-stable RK method. 
Hence, if the RK method 7Z S is A-stable, in order to preserve the hyperbolic 
structure, we only need to require 

Now, we show that the above inequality can be equivalently reduced to 

(2.25) det(J - zA) > 0. 

Indeed, since 7Z S is A-stable, the eigenvalues of its coefficient matrix must 
lie on the right half plane, which implies det(/ + zA) > for z > 0, and this 
together with (|2.24p further implies that det(J — zA) > 0. Next, using again 
the fact that the eigenvalues of A lie on the right half plane, we have that 

det(J + zA) 
det(/ - zA) > ' 

if z > satisfying det(I — zA) > 0. Therefore, the hyperbolic equilibria 
structure-preservation region of an A-stable SRK method TZ S is given by 

(2.26) {h > : det(J - zA) > 0, z = f3h}. 

For example, the hyperbolic structure-preservation region for the well-known 
midpoint formula, where A = 1/2 is given by 

(2.27) {h > : 1 - -f3h > 0} = {h > : 0h < 2}. 

It is readily seen that even for SRK method which possesses good classical 
stability properties, we should still require some restrictions on its step-sizes 
for practical computations. Therefore, the investigations on the equilibria 
structure preservation provide a novel and useful criteria to choose step- 
sizes for symplectic integrators in addition to the classical linear stability 
requirements. 

We now apply an s-stage SPRK method TZg — TZg to system f)l . 16j) and 
get the following scheme 

_ flP) (-*)-! 



(2.28) 




h PJi 

h<Ui 
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where R®(z) = det(7 - zAW + ze s b^ T )/ det(J - *AM)(i = 1,2), are the 
stability matrices for lZg(i = 1, 2). Consequently, the hyperbolic structure- 
preservation region for TZ S — is 

(2.29) {h>0: R {1 \-z) < 1 and R {2) {z) >l,z = ph}. 

For the well-known Lobatto IIIA-Lobatto IIIB pair, since both 

Lobatto IIIA method (TZg ) and Lobatto IIIB method (7^s ) are symmet- 
ric A-stable RK method, and therefore R®(z) = det(J + zA (i) )/det(I - 
zA®)(i = 1,2). Similar to (I2T26I) . condition (pT29j) is reduced to 

(2.30) {/i > : det(J - zA {i) ) > 0, i = 1,2, z = 

for Lobatto IIIA-Lobatto IIIB method, where in particular, we note that 
det(J - zA®) = det(J - zX AW ) = det(J - zX A(2) )(i = 1,2) with X AW (i = 
1, 2) given in (12TT61) and OTH) . 
In summary, we have 

Proposition 2.2. T/ie hyperbolic equilibria structure-preservation region of 
an SRK method 1Z S for system U.16\) is 

det(I-zA) , det(I + zA) 

{h>0: TT < 1 and ^-77 tt >1,Z = /%}, 

1 det(/ + zA) det(J-^A) ' K J ' 

and i/izs condition is reduced to 

{h>0: det(J - zA) > 0, z = /%}, 

z/7^ s is A-stable. For an SPRK method, 1Z S — TZa , the hyperbolic structure- 
preservation region is given by 

{h>0: R {1 \-z) < 1 and i? (2) (z) >l,z = 0h}, 
where R^ and R^ are, respectively, the stability functions for and 

Now, we consider the structure-preservation of the symmetric composition 
method (|2.20p . The 4th order symmetric composition of mid-point formula 
for system (|1 . 16[) is 



(2.31) 



(1 - 7iz)(l - 7 2 z)(l - 73*0 
" f l +7iz)(l +7 2 2;)(1 +732) 



1 +7iz)(l + 72*0 1 1 + 73*) 
= n VTi ^Fi \ « 



(1 - 7iz)(l - 72*0(1 - 73*) 

with z = (3h, 7i = 1/(2 - 2 1 / 3 ) and 72 = 2 1 / 3 /(2 - 2 1 / 3 ). Therefore, for 
scheme to preserve the hyperbolic structure, we have 
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which yields 

2 - 2 x /3 

(2.32) < z < (sa 0.587401). 

Clearly, in comparison with (I2,27h . the structure-preservation region be- 
comes smaller after composing. This result applies equally to the SPRK 
method. Consider the following symplectic Euler method for system (|1.16p . 

{Pi+i = (1 - z)Ph 
1 RU 
m+i = Y^~z qu z = 

whose structure-preservation region can be verified to be 
(2.34) < z < 1. 

Since the 2nd order symmetric composition of symplectic Euler method is 
the midpoint formula, hence the structure-preservation region of the 4th 
order symmetric composition of scheme (12.33H is given by (|2.32p . which is 
obviously decreased in comparison with (j2.34[) . 

3. Extensions to nonlinear Hamiltonian systems and some 

applications 

In this section, we extend the analyses in the previous sections to the 
nonlinear 1-d.o.f. Hamiltonian system 



(3.1) 



P = -dH/dq := f(p,q), 
q = dH/dp := g(p,q), 



with (p, q) G R 2 . We denote by $n the set of equilibrium points for sys- 
tem (|3.ip . i.e., 

Sh ■= Up,q) G M 2 ; V M H(p,q) = oj. 

The stability matrix for system (|3,ip is locally defined for every (p, q) E &h 
as 

df/dp df/dq 
dg I dp dgj dq 



(3-2) 3s(p,q) = 

whose eigenvalues are given by 



(p, q), 



Ai, 2 = ±VD 

with 

D:=D{P > q) = [ Tqdp-dp-dq- ]{P > q) - 
Therefore, we can make the following classification: 

(i) If D(p,q) < for all (p, q) E (oh, then the equilibrium points for 
system (13.1j) are all of elliptic type, and we then call system (|3.ip 
having elliptic equilibria structure; 
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(ii) If D(ft,q) > for all (p,q) E $h, then the equilibrium points for 
system (13. If) are all of hyperbolic type, and we call system (|3.ip 
having hyperbolic equilibria structure; 

(hi) If D(p,q) > for some (p, q) E S'h, while D(p,q) < for some 
(p,q) E $h, we call system (13. ip having mixed- type equilibria struc- 
ture. 

As an example, we suppose that system (I3.ip have hyperbolic structure 
and apply an s-stage SRK method 1Z S for its discretization to yield 

f P = Pl e s + hAF(P,Q), Pl+1 = Pl + hb T F(P,Q), 
I Q = qie s + hAG(P, Q) , q l+l = q t + /i& T G(P, Q) , 



(3.3) 



where F(P,Q) = [f(P l ,Q 1 ),...,f(P s ,Q s )) T and G(P, Q) = [g(P 1 ,Q 1 ), . . . , 
g(P s ,Q s )] T are the internal stage vectors. The corresponding discrete dy- 
namical system for scheme (13. 3p is 



(3.4) 



5+ Pl = b T F(P,Q), 
5+ qi = b T G(P,Q). 



'N 



Clearly, the points in Sh are still the equilibrium points (or fixed points) of 
(ET3j) . The stability matrix for (pT4"|) is 

b T dF/dpi b T dF/dqi 
b T dG/dpi b T 3G/dqi 

Therefore, in order to preserve the underlying equilibria structure for the 
numerical scheme, we need to study the matrix structure of 3n(P,q) with 
(p,q) E &H-> as what we have done before. Obviously, such arguments apply 
equally to SPRK methods, and to Hamiltonian systems having the other 
two kinds of equilibria structures. In the sequel, as an application, we give 
an example. Consider the following separable 1-d.o.f. Hamiltonian system 



(3.5) 



P 



ap(l 
a(2p 



P), 
1)9, 



where a > 0. This problem has been studied in [3], and the solutions p(t) j 
l(i — > +oo) if < p(t ) < 1, and p(t) [ l(t —> +oo) if p(t ) > 1; while for 
any q(t ) > 0, q(t) ) +oo(t -> +oo) if p(t ) > 1/2, and if < p < 1/2, q(t) 
is first monotonically decreasing for t < to = (ln(l — p(to)) — lnp(to)/ Ina) 
and then q(t) | +oo for t(> to) +oo. It can be seen that (0,0) and (1,0) 
are two hyperbolic equilibrium points for (13. 5p . The stability matrices here 



arc 



(3.6) J S (0,0) 



O: 








-O: 



and J S (1,0) 



-a 
a 



The symplectic Euler method for (|3.5p is read as 

j pi +x = (1 + z)pi - zpf, 
) Ql+l = Ql + z(2pi - l)q i+1 , 



ah. 
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or, equivalently, 

dfpi = api{l-pi) 
(3-7) { x+ (2p,-l) 



1 + z — 2zpi 

We have the following theorem for scheme (|3.T|) to give a true simulation 
(see m) 

Theorem 3.1. The sequence {p{\ — > 1 and {qi} — > +oo osn-> +oo i/jf 

1 + z 

(3.8) < z < 1 and p(t ) < . 

2^ 

If condition (|3.8p is destroyed, then it is shown in [3] that spurious solu- 
tions or even periodic solutions will be encountered. Through straightfor- 
ward calculations, the stability matrices of the discrete dynamic (13.70 at its 
equilibrium points are given by 



Jjv(0,0) 



a 

-«/(! + z) 



and Jzv(l,0) 



-a 
a/(l-z) 



In comparison with the stability matrices in (13. 6p , we see that for scheme ([37 
to preserve the hyperbolic equilibria structure, we need require 1 — z > 0, 
i.e., < z < 1. That is, the preservation of equilibria structures only endows 
the numerical integrator a prerequisite for successful simulations. 

For Hamiltonian systems of higher dimensions, the situation will become 
much more complicated. However, we can make use of the usual linearization 
techniques for our investigations, which are interesting topics for our future 
study. 
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